/*
 *  MWSolve.cpp
 *  T3nsors
 *
 *  Created by Michael Barriault on 10-07-12.
 *  Copyright 2010 University of Guelph. All rights reserved.
 *
 */

#include "MW.h"
#include <cmath>
#include <sys/time.h>
#include <iostream>
#include <fstream>

void MW::Solve(double tn, double s, int icid) {
	std::ofstream info((id + "/info.txt").c_str());
	info << id << std::endl;
	info << "n\t" << O[0]-1 << std::endl;
	info << "e\t" << epsilon << std::endl;
	info << "k\t" << K << std::endl;
	info << "c\t" << c << std::endl;
	info << "w\t" << w << std::endl;
	info << "s\t" << s << std::endl;
	info << "T\t" << tn << std::endl;
	info << "j\t" << J << std::endl;
	info << "damp\t" << damp << std::endl;
	info << "p\t";
	FORa info << O.p(a) << " ";
	info << std::endl;

	Badge->PrintStatus("Forming solution.");
	double dt = s*O.d(0);
	int N = int(tn/dt);
	double t=0;
	double I = J;
	int z0 = (O[2]-1)/4;
	
	int MemSize = 0;
	MemSize += 3+2*3; // x/Phi/Psi
	MemSize += 10+10*3; // RK
	MemSize += 9; // dots
	MemSize += 2; // R,dR
	MemSize += 5; // Variables sent to Print
	MemSize *= 64*O();
	
	Badge->PrintStatus((string)"Estimated memory usage: " + bytes(MemSize));
	
	long long int StoSize = 13LL*8*O()*N;
	
	Badge->PrintStatus((string)"Estimated storage usage: " + bytes(StoSize));

	Scalar R(O), dR(O);
	FOR(i,O[0]) FOR(j,O[1]) FRO(k,z0,3*z0) {
		R(i,j,k) = smoother(x(2)(i,j,k),1);
		dR(i,j,k) = 8*x(2)(i,j,k)*pow(pow(x(2)(i,j,k),2)-1,3);
	}
	
	Badge->PrintStatus("Initial data");
	if ( icid >= 0 and icid < 10 ) { // plane-wave
		FORO {
			E(0)[o] = R[o]*cos(w*x(2)[o]);
			E(1)[o] = R[o]*sin(w*x(2)[o]);
			B(0)[o] = dR[o]*cos(w*x(2)[o])/w - R[o]*sin(w*x(2)[o]);
			B(1)[o] = dR[o]*sin(w*x(2)[o])/w + R[o]*cos(w*x(2)[o]);
			if ( icid == 1 ) E(2)[o] = R[o]*c; // divE violation
			if ( icid == 2 ) B(2)[o] = R[o]*c; // divB violation
		}
		info << "E0 = R cos(wz) e_x + R sin(wz) e_y";
		if ( icid == 1 ) info << "+ R c e_z";
		info << std::endl;
		info << "B0 = ( dR/w cos(wz) - R sin(wz) ) e_x + ( dR/w sin(wz) - R cos(wz) ) e_y";
		if ( icid == 2 ) info << "+ R z e_z";
		info << std::endl;
	}
	if ( icid >= 10 and icid < 20 ) { // simple polynomial
		FORO {
			E(0)[o] = R[o]*smoother(x(1)[o],0.2);
			B(1)[o] = R[o]*smoother(x(0)[o],0.2);
			if ( icid == 12 ) E(1)[o] = R[o]*c; // E.B violation
		}
	}
	if ( icid >= 20 and icid < 30 ) { // all-sine
		FORO {
			FORa E(a)[o] = R[o]*sin(w*x(2)[o]);
			B(0)[o] = E(0)[o];
			B(1)[o] = -E(0)[o];
		}
	}
	info.close();
	
	Badge->PrintStatus("Evolution loop");
	struct timeval start,first;
	struct timezone tzp;
	gettimeofday(&start,&tzp);
	for ( int n=0; t<tn+dt/2 and !stopit; n++ ) {
		if ( fabs(t-floor(t)) < dt )
			J = I/(l2(E*E+B*B)+1);
		gettimeofday(&first,&tzp);
		if (all) {
			Print(t, Phi);
			Print(t, Psi);
			Print(t, E);
			Print(t, B);
			Print(t, E*B);
			Print(t, B*B);
			Print(t, div(E));
			Print(t, div(B));
		}
		Print(t, E*E);
		Badge->PrintTime(t);
		t += dt;
		RK3(dt);
		Badge->PrintStatus((string)"Time remaining: " + timer(first,N-n),false);
		Badge->Progress(t/tn);
	}
	Badge->PrintStatus((string)"\nTotal elapsed time: " + timer(start));
}

